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Abstract 

We investigate an operator renormalization group method to ex- 
tract and describe the relevant degrees of freedom in the evolution of 
partial differential equations. The proposed renormalization group 

approach is formulated as an analytical method providing the fun- 
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damental concepts of a numerical algorithm applicable to various 
dynamical systems. We examine dynamical scaling characteristics 
in the short-time and the long-time evolution regime providing only 
a reduced number of degrees of freedom to the evolution process. 
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64.60.Ht 

KEY WORDS: Renormalization group; coarse-graining; nonlinear evo- 
lutionary dynamics; partial differential equations. 

1 Introduction 

Based on the arguments of L.P. Kadanoff jy, K.G. Wilson developed a 
(numerical) renormalization group (RG) approach in the context of critical 
phenomena [|| and later used the developed formalism to solve the Kondo 
problem ||. fn practical RG applications to systems with many degrees 
of freedom a RG transformation is established for a stepwise elimination 
of the irrelevant degrees of freedom in the system. The universality of the 
critical behaviour of the physical system emerges in the RG approach in a 
natural way, with a set of universal exponents for each universality class. 
Soon after the considerable success of the RG for equilibrium systems a 
generalization of the method was introduced to handle dynamic properties 
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of spin systems Q. The goal was to obtain critical exponents for stochas- 
tic equations and the developed generalization of the Wilson's equilibrium 
theory has become known as the dynamic RG (DRG) method for non- 
equilibrium systems [BJ . The DRG method, inherently defined as a Fourier 
space technique, was further developed and improved by other authors in- 
cluding D. Forster, D. Nelson and M. Stephen who studied the dynamics of 
the noisy Burgers equation In 1986 Burgers equation was transformed 
into the language of growth processes by M. Kardar, G. Parisi and Y.-C. 
Zhang later called the KPZ equation. Applying the DRG to the KPZ 
equation the corresponding universal critical exponents were derived per- 
turbatively within various higher order expansions ^ ^| . 
However, the impact of dimensionality on applying the DRG method to 
the KPZ or the related Burgers equation is crucial JTol. For dimension 
d > 1 the exhibited strong coupling RG flow behaviour is not accessible 
by perturbative expansions and the standard approaches fail to calculate 
critical exponents pj]] . New analytic approximations were proposed |]l2f 
including RG methods based on technical variations or combining different 
approaches like perturbation and scaling techniques []l3| . Furthermore nu- 
merical simulation techniques were developed to allow for estimating the 
critical exponents |l4|, [l5| although these approaches can suffer from nu- 
merical instabilities [jl6| . 

Contrary non-perturbative Real-Space RG techniques have been proven to 



produce accurate results within strong-coupling regimes ]17, 18 and there- 
fore offer, in principle, the possibility to study models that belong to the 
KPZ universality class. Why is it then that currently no satisfying gener- 
ally applicable Real- Space RG approach is available for dynamical systems 
as this is the case for their equilibrium counterparts? In the following we 
would like to answer this question in detail and we propose a list of the 
fundamental problems arising in the context. 

Many of the modern Real-Space RG (RSRG) methods for strong-coupled 
systems use an operator formalism to define the RG transformation (RGT) 
[ fi9[ pc| which under iterative application defines the RG flow. The RGT 
in the operator formalism is formally based on the application of two linear 
maps, a truncation or fine-to-coarse operator and an embedding or coarse- 
to-fine operator |2lJ . Usage of these two operators should be sufficient 
for the construction of renormalized or effective quantities. Within the 
equilibrium operator RSRG an effective observable, including in particular 
the Hamiltonian, is defined by first applying the coarse-to-fine operator 
followed by the original observable. Finally the fine-to-coarse operator is 
used as a map back on the effective vector-space. 

Recently N. Goldcnfcld, A. McKane and Q. Hou investigated the use of 
operator RSRG methods to solve partial differential equations (PDEs) nu- 
merically p2| and carried on their ideas in a further paper pjJ. In their 
work the authors generalized the operator concept of the RSRG to non- 



equilibrium system by replacing the equilibrium Hamiltonian operator by 
the time evolution operator of a partial differential equation. However the 
authors showed that such an approach is, in principle, not well-defined. 
Furthermore the proposed approach is less systematic and unsatisfactory 
for several reasons. Here we summarize the concerns of Goldcnfeld et al 
and complete the list by considerations of the present authors. 

1. Neither the coarse-to-fine operator nor the fine-to-coarse operator 
are uniquely defined. The chosen coarse-graining procedure does not 
depend on the specific dynamical system. 

2. The coarse-to-fine operator and the time evolution operator do not 
commute in general. Therefore no consistent and well defined RSRG 
operator formalism can be established. 

3. The geometric construction of the coarse-to-fine operator is unsatis- 
factory since it assumes the relevant degrees of freedom to be dis- 
tributed within the long wavelength fluctuations. This is not neces- 
sarily the case for the evolution of many degrees of freedom interact- 
ing through nonlinear dynamics. 

4. Instead of coarse-graining the governing PDE in a systematic proce- 
dure, the provided initial field configuration at time t = is coarse- 
grained with respect to a larger lattice spacing. Evolving a coarse- 
grained field configuration is not equivalent to successive RGTs re- 
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moving the relevant degrees of freedom in time. 

5. In a practical application of the operator concept introduced by Gol- 
denfeld et al the coarse-graining procedure and the time evolution 
procedure are successive operations. However, performing an equi- 
librium RG step followed by an evolution of the system in time is, in 
general, not equivalent to a dynamic RG step. 

6. The proposed operator scheme does not allow to calculate the dis- 
tribution of the relevant degrees of freedom in a dynamical system. 
Accordingly, RG flow equations do not yield any insight to the rele- 
vant physics of the system. 

In our previous work pjfl , which can be considered as a further develop- 
ment of the work of Goldenfeld et al J23], we developed a consistent and 
mathematically well defined RSRG operator approach. The basic concepts 
also used in this article are briefly summarized in section ^ and general- 
ized by including time evolution operators. In our previous work we used 
the introduced operator formalism to construct the embedding and trun- 
cation maps by a quasi-static geometric coarse-graining. The conceptional 
disadvantage of a geometric approach is that it inherently integrates out 
the small scales which, because of scale interference in the evolutionary 



process, are not necessarily the scales one wishes to ignore 22 . However, 



the technique can be easily applied to a great variety of physical problems, 



even if defined on very large lattices. 

The essential difference in the approach proposed here is a non-geometric 
construction of the coarse-to-fine and fine-to-coarse operators involving the 
evolution dynamics of the particular PDE of interest. In this work we solve 
all reported problems summarized above and calculate observables to de- 
termine the universal characteristics in a class of nonlinear PDEs. As a 
non-geometric generalization of our earlier RSRG operator approach the 
proposed method provides a non-perturbative real-space RG analogue to 
the DRG Fourier space method. 

We continue to present the fundamental concepts of the non-geometric 
RSRG operator approach in section ^|. In particular, we develop two con- 
cepts applicable in the short and long time regime of the evolutionary 
dynamics respectively. In section ^ we solve the linear diffusion problem 
exactly and present an analytic construction scheme for the coarse-to-fine 
and fine-to-coarse operator. We continue in section |^ by generalizing the 
method as a numerical algorithm applicable to nonlinear evolution equa- 
tions in general. Using transformations which are local in time allows for 
exact mapping of the non-linear dynamics to linear evolution dynamics 
governed by the diffusion equation. 

In this work we denote the dependence of a function / on spatial variables 
x and temporal variables t as f(x,t). Contrarily the spatial and temporal 
dependence of a discretized function / will be denoted by indexing sets 
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i, i + 1, . . . and t, t + 1, . . . respectively. Here i, i — 1 and i,i + 1 refer to 
neighbouring sites in the lattice whereas t — l,t,t+l denote successive time 
steps. Not necessarily neighbouring sites are denoted as i,j and by using 
capital letters /, J we refer to the sites of the effective lattice. Quantities 
defined in this effective vector space are equipped with a prime. A set of 
vectors is denoted as {|vi)}fl M where the indexing set starts at i = M and 
extends to i = N. 

2 The Operator Real-space RG Approach 

In this work we consider examples of evolution equations of the form 



with <f> — t) a function of space x and time t and the operator H acting 
as the generator of the evolution. Although every operator H acting locally 
in time can be considered, we restrict our calculations for clarity to linear 
and quadratic evolution operators, H and Q respectively, 



d t <f> 



Hcj>, 



(1) 



H 



(H + Q). 



(2) 



Discretizing equation (|l|) in space and time using (|^) yields 



N N 



<t>i,t+At = + Af • H it j)<pj >t + At Qi,j,k<Pj,t<t>k,t •= fi[<Pt] 




(3) 
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with i £ {1, . . . ,N} and N denoting the number of sites in the lattice. 
The spatial lattice spacing is denoted as Ace and the discrete temporal 
integration interval as At. The function <f> t defines a vector in a vector 
space V N and the functional / [<f>t] defined in (|^) is a map 

/ [<j> t ] : V N C V N . (4) 

Using this vector space notation the concept of a truncation or fine-to- 
coarse operator G together with an embedding or coarse-to-fine operator 
G p is provided by the not necessarily commuting diagram 

AT r AT /r\ 

V N c V N (5) 

G g p and M < N . 

Here the truncation operator G is defined by the relation 

<t>'i = y] Gi,t & , (6) 

i=l 

where capital indexing letters / £ {1, . . . , M} refer to lattice sites in the 
effective vector space V M . If M = N and Ker(G) = 0, the embedding 
operator G p may be chosen as the natural inverse of the truncation operator 
G defined as 

M M 

4>i = = E G a^ ( 7 ) 

i=i i=i 

and the diagram (^) commutes. 

However, within practical applications the idea is to choose M < N and 



replace equation (g) with a coarse grained evolution equation providing less 
degrees of freedom. Considering G as a projector on the relevant degrees of 
freedom the inverse operator G _1 does not exist. This naturally demands 
for a generalized definition of the embedding operator G p as the pseudo- 
inverse 

M 

& = £G^'z (8) 
1=1 

calculated by Singular Value Decomposition (SVD) p5| . For M < N the 
functional /' [<j)' t ] is defined as in (||) on an effective lattice with a reduced 
number of lattice sites M 

f [<j>' t ] : V M C V M 

M 

and (j>' I<t I <t>'i, t+ At = i 11 ^ + At ■ H 'i,j) 

M 

+ At ]T Q'j,j, K (t>'j^'K,t ■ (9) 

J,K=1 

In equation (||) H' and Q' denote the effective linear and quadratic evolu- 
tion operators defined as Q 

:= GzjHijG^j and Qi,j,k '■= Gi,iQi,i,kGjjG kK . 

(10) 

Inserting (||), (|J) and (0) into the non-commuting diagram (||) we replace 
equation (|3j> by an approximate field evolution equation on a coarse grained 
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lattice as 

N M 

<Pi,t+At « J]E (! ?,l( 1 V + At, W^ 
i=l ,7=1 

JV M 

+ X! G llQ'l,J,K G J,j G K,k(t>j,t<t>k,t 

j,k=l J,K=1 

G p (l + At-H' + At-Q')G(t> t ] , (11) 

- i 

describing the evolution of the field <fi under a reduced number of degrees 
of freedom. 

Equation (|ll|) defines a RG transformation (RGT) within the operator 
formalism and iterating the RGT defines a RG flow. Carrying out one 
RGT is called a RG step (RGS) and according to the concept developed 
in this section is equivalent to an approximate field evolution using less 
degrees of freedom. Equation (|ll|) defines the real-space analogue of a 
RGT established within the DRG method. The field itself provides the set 
of parameters used to establish a RGT pp| . Furthermore, equation ( p"T| ) 
fuses the coarse-graining and the time evolution procedure and provides a 
solution to the problem number 5 formulated in the introduction. 
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3 Non-geometric reduction of the degrees of 
freedom 

In this section we provide a general concept for reducing the degrees of 
freedom in evolutionary systems without geometrically coarse-graining the 
lattice equations. Including the temporal characteristics of the particular 
partial differential equation (PDE) into the construction of the embedding 
and truncation operators we distinguish between a short-time regime and 
a long-time regime. 

3.1 Short-time evolution 

Within the short-time regime we describe the evolution of the field as a 
perturbation of the initial field configuration at time t = 0. Evolving the 
initial field </>o for M time steps At the dynamics within this short-time 
interval is conserved by the set of vectors 

S := {0o, H<f>o,H 2 <p , H M - X M (12) 

where H is the evolution operator defined in (Q). Using the set of vectors in 
( [l2| ) as the columns and rows of the linear operators G and G p (previously 
orthonormalized) , these can be considered as projection maps from and 
into the space V M respectively. 
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If M < TV the relation 



\ M 

ih = G p il + AtH') G<po with H' = H' + Q' (13) 



governs an exact evolution of the field <j>t on the effective coarse-grained 
lattice for t < MAt. In this case the states in (|l^ ) span a subspace V M 
of the full vector-space V N conserving the relevant degrees of freedom for 
the short-time evolutionary regime. In the considered short-time regime, 
we may rewrite equation ([l^) as 

r / - \ i M 

<t> t = G p \1 + AtH'J G (j} (14) 

which determines a RG flow in the short time regime and a RGS is defined 
by relation (|Tl|). 

However, if t > MAt equation (|l4|) is an approximation to the evolved 
field 4>t- In this case the approach is only applicable if no relevant scale 
interference in the evolutionary process occurs. This is unlikely the case 
for longer times in nonlinear dynamical systems and relation (|lj) becomes 
a crude approximation giving rise to numerical instabilities. 



3.2 Long-time evolution 

Nonlinear evolution processes exhibit most of their characteristics in the 
long-time regime. The asymptotic form of the field or surface configuration 
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in growth phenomena H or the formation of turbulent states out of spiral 
waves p4] are only two examples. This gives rise to a reduction scheme 
for the degrees of freedom of a dynamical system valid in the long-time 
asymptotics. This in turn demands for new concepts in the construction of 
the embedding and truncation operators G p and G away from any initial 
field configuration. 

Therefore our task is to minimize the magnitude of the non-commutativity 
in the diagram (jE)) for all times of the evolution process. This in turn 
rises the important question if such a minimum exists and if a numerical 
algorithm will uniquely converge. To measure the non-commutativity in the 
diagram (^) we calculate the difference between the approximately evolved 
field using a reduced number of degrees of freedom and the exactly evolved 
field as 



We call £ the error operator and generalize the operator concept introduced 
in section || according to the commutative diagram 



£<j> t := (l + AtH)<p t - (G p )(l + MH')G4> t . 



(15) 




■N 



and M < N 



(16) 
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The introduced minimization concept is one of the fundamental principles 
in equilibrium operator RSRG techniques |l7], pT| . Using the established 
notation in the formulation of the dynamic operator RSRG, (1 + AtH)<f>t 
is called the target state jlTj and (G p )(l + AtH')G <pt an optimal repre- 
sentation of the target state jl7| in the subspace V M C V N . 
We are interested in the limit £(j) t — > for t 3> subject to any field config- 
uration 0o- To exclude the explicit field dependence from the minimization 
procedure we rewrite equation ( filf ) in operator form as 

£ = (1 + Atff) - G p (l + AtH')G . (17) 

To make this operator as small as possible we minimize £ in the matrix 
notation according to the Frobenius norm [|]. Inserting the definitions ( [To| ) 
we rewrite equation ( |l7| ) using the Frobenius norm \ \p as 

\£\ F := {l-(G p G) + At\ i H-(G p G)H(G p G)]y . (18) 

The measured magnitude of the non-commutativity in the original diagram 
(0) has become a number \£\p greater or equal to zero. According to 
equation (|l8j ) the minimization of the error operator £ only depends on 
the composed operator G P G. Therefore two different G matrices for which 
G P G is the same operator V N c >■ V N yield the same error number 



tin principle every matrix norm can be used, although the Frobenius norm can be 
easily related to concepts from linear algebra like singular value decomposition. It is 
defined to be the sum of the squares of all the entries in a matrix. 
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The matrix G P G is called the reduction operator since it governs the field 
evolution under a reduced number of degrees of freedom. To make this more 
obvious we rewrite G and G p in terms of a Singular Value Decomposition 

M M 

G = ^AiWH and G p = £ A^ 1 |v 4 ) (u. ( | , (19) 

i=l i=l 

where {Ivj)}^ and {luj)}^ are respectively sets of iV-dimcnsional and 
M-dimensional vectors. By means of ( p"9"| ) wc define the reduction operator 
as 

M 

n := G P G = X>i><Vi| • (20) 

i=l 

In this notation the reduction operator 1Z is composed of M states and acts 
on the original vector space V N , with N > M. Since there are some states 
{l v i)}i=M+i which together with {x~i}f£ 1 represent a whole orthonormal 
basis of Vmi the reduction operator 1Z projects out those states representing 
the less relevant degrees of freedom in the evolution process. 
We would like to point out that different truncation operators G can be used 
to construct the same reduction operator 1Z = G P G. In fact every RSRG 
approach requires a correct choice of the RG transformation (RGT) which 
is not uniquely defined. In the operator formalism the embedding and 
truncation operators are used to construct the RGT p(J. In the dynamic 
operator RSRG method introduced in this work, the RGT is defined as 
one time step in the field evolution equation (|ll]) . From this point of view 

1G 



all truncation operators yielding the same reduction operator 1Z represent 
an error equivalence class as far as the evolution error is the only concern. 
Other computational considerations, such as the sparseness of the H matrix 
shall force us to choose one or another truncation operator. 
Before we discuss the field evolution according to the reduction operator 
1Z we would like to give a remark concerning the structure of equation 
Choosing At = 0, i.e. no dynamics is involved, equation ([t8]) reduces 
to the geometric operator real space RG approach pij. For At > the 
third summand on the right hand side includes the evolution operator of 
the particular PDE into the minimization process. Furthermore this term 
is weighted by the time interval At. Therefore the particular evolution 
operator is included into the minimization process and resulting relevant 
degrees of freedom for the evolutionary dynamics can depend on the chosen 
time interval. 

To reformulate the approximate field evolution in terms of the reduction 
operator we insert the definitions ([To]) and (|20| ) into (|ll|) resulting in 

&,t+At = n[4> t + At-{H + Q)n^ t ] . (21) 

Within practical calculations the process of constructing 1Z can be decom- 
posed according to successive dimensional reduction of the target space 

yN c ^ yN-1 c ^ yN-2 c ^ , , , C ^ yM _ ( 2 2) 
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Within the first dimensional reduction V N c >- V N ~ X we have to mini- 
mize ([l8|) by adjusting the N components of the vector |vjy) which defines 
the first order reduction operator as 

N-l 

II 1 := l v *> ( v *l = 1 - Wn) (vjvI , (23) 

i=l 

where we used the notation of relation (po|) . Analogously we calculate the 
second order reduction operator as 

N-2 

n 2 := ^l v *X v *l = i - l v w) ( v ^vl - l v w-i> ( v w-il 

i=l 

= K 1 -V\ VN _, ) . (24) 

In ( p4] ) we have introduced the notation of the projection operator 
defined as the projection map V\ v \ = (|v) (v|) onto the state |v). 
Iterating up to the Mth order the reduction operator is calculated as 

u m = n M+i _ Wn+1 _ m) ( Vjv+1 _ m | = n M+1 -V\ VN+1 _ M) , (25) 

with 1Z° — 1. In equation ( p5| ) we have denoted by V\ VN+1 _ M \ the projec- 
tion operator onto the state |vjv+i-m) which is the target to calculate in 
the Mth minimization procedure. 

Similar calculations can be performed for the error operator £. Using the 
abbreviation |v) := |v)jy the error operator for the first dimensional reduc- 
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tion is defined as 



E 1 = l-(l-\v){v\)+M H- (1 - \v) (v\)H (1 - \v) (v\ 



V\ v) + At(HV\ v) +V\ V) H- (v\ H \v) V\ v) 



(26) 



Using the iterative scheme ( p5| ) we write the minimization procedure for 
the state |vat + i_m) as 



\£ M \ F = 1-K M + At( H - K M HTZ 



(27) 



According to relation (27) we call a truncation V N c y N ~ M of the 

original vector space V N an order M truncation or a reduction of order M 
in the degrees of freedom. The ratio 



A := N/(N-M) 



(28) 



is defined as the reduction factor A. 



4 An Exactly Solvable Lattice Model 

In this section we treat the special case H = H, i.e. a linear evolution 
dynamics. For both a selfadjoint and a non-selfadjoint evolution operator 
the reduction operator can be calculated exactly up to an arbitrary order in 
dimensional reduction of the original vector space V N . As an example for 
a selfadjoint linear evolution operator we examine the diffusion equation 
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given by 



dw(x, t) 
dt 



v ■ V 2 w(x,i) , 



(29) 



where the diffusion coefficient v describes the strength of the relaxation 
process. 

In figure [l] we evolved the initial field configuration (dashed line) according 



to the method described in section 3.1 for a total time t = 100 and t — 1000 



The solid curve shows the field evolution on the original lattice composed of 



Field evolved with ell the degrees of freedom 

Field evolved with half of the degrees of freedom □ 
Initial field configuration — 




(b) 



Field evolved with all the decrees of freedom 

Field evolved with, half of the degrees of freedom □ 
Initial field configuration ---- 




20 25 



Lattice site s 



Figure 1: The initial and the evolved fields for a lattice composed of 32 
sites after 100 and 1000 time steps (a) and (b) respectively. The initial field 
configuration is drawn by the dashed curve. The fields were evolved using 
a time interval At = 1 and a diffusion parameter v — 0.01. The numerical 
minimization was performed up to the 16th order reduction in the degrees 
of freedom equivalent to a reduction factor of A = 2. 

32 sites. The field evolved with half of the degrees of freedom corresponding 
to an evolution on a lattice composed of 16 sites is displayed by a □ for 
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every lattice site. The embedding and truncation operators G p and G were 
composed according to ( p"2| ) with M=16. For the numerical evolution of 
equation (|29|) we consider the following discrete differencing of (J2S|) 



At 



Wj+i,t - gw^t - w. t _ M 
(At) 2 



(30) 



In equation (30) we used a forward Euler scheme in time |25| and the linear 
evolution operator H is defined as 



:= < 



-2u/(At) 2 if i = j 



f / (At) if j nearest neighbour site of i 



else . 



(31) 

From figure [l] it is obvious that the evolution under a reduced number of 
degrees of freedom does not differ from the evolution of the field using all 
degrees of freedom. Even for longer evolution times which corresponds to 
equation ( |l4| ) with M ^> N the field evolved with half of the degrees of 
freedom does not differ from the exact evolution. This indicates that the 
relevant degrees of freedom can already be deduced from the early time 
dynamics. Therefore the degrees of freedom do not mix in linear evolution, 
i.e. small scale dynamics and large scale dynamics do not interfere. 
To calculate the long-time evolution characteristics we have to minimize 
\£\f in ( |l8| ) by iterative application of ( |27| ) according to (p5|). We di- 
agonalize the self-adjoint operator H according to the orthonormal basis 
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{| M i)}i=i with real eigenvalues Ei given by 



JV 



(32) 



i=l 



The eigenvalues Ei are supposed to be real and in increasing order: E\ < 
Ei < . . . < En. If the basis {|uj)}£Lj is orthonormal the vector \v) in ( p6| ) 
can be decomposed according to 



3 = 1 



(33) 



where we are assuming fij £ R for all i. Using relation (53) the different 



summands in equation (|26j) can be calculated explicitly as 

TV N 

P\v) = ViVj \ui){uj\ (v\H\v) = Y^ El = W 



N 



i=l 
N 



H ^\v) = Y ^0 E i \ U i) 



and relation (E6J) becomes 



£ij = t + At (Ei + Ej - (H)) fXifij 



(34) 



(35) 



According to definition (18) we have to minimize 

N 



]T [l + At(E t + E ] -(H)) 

i,j=l 



1 1 
ViVi 



(36) 



Expanding into explicit summands and retaining only terms which are first 
order in At this yields 



\£\ F = 1 + 2At(H) 



(37) 
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Therefore (H) must be minimized as a function of the restricted to 

the condition ^ /i| = 1. To include the constraint in the minimization we 
introduce a Lagrangian parameter A as 

N I N \ 

(H) = ^-1 • (38) 

»=1 \i=l / 

Derivation with respect to each of the \ii yields 

(Ei + X)fj,i = for all i e {1, . . . , N} . (39) 

Thus, for each i, either /ij = or Ei = A. But, if all eigenvalues Ei are 
different the latter condition may apply to only one value of i. In order to 
fulfill the unit norm condition, the i-th component must have value /i, = 1 . 
Therefore the state \v) which minimizes the error operator is 

/ii = 1 and [M = if % ^ 1 , (40) 

because it corresponds to the smallest eigenvalue of H . In the case of 
degenerate eigenvalues, there is not a unique minimum vector, but a mini- 
mum eigen-subspace. 

To validate this result numerically we minimized the error operator ( p7|) 
of order M = 8 and M = 16 by calculating the set of states {|vi)}^ =Ar _ M 
iteratively. To perform the minimization we used a sequential quadratic 
programming (SQP) method |^6| provided by the NAG Library |27]]. The 
SQP algorithm used the 32 components of each of the vectors in the set 
{l v i)}^7V-M as adjusted variables and the necessary 32 gradient elements 
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have been numerically estimated within the NAG Library. Both for a 
truncation of order M = 8 and M = 16 the SQP algorithm converged 
into a unique minimum stated by the NAG Library as the optimal solu- 
tion found. The computer time required to evaluate a reduction operator 
of order M = 16 was about an hour on a SUN-SPARK Ultra 10 work- 
station (SUN Microsystems, Inc., Palo Alto, USA). In agreement with the 
analytical scheme introduced above we always recalculated the eigenstate 
corresponding to the highest eigenvalue within machine precision. 
In figure || we plotted the corresponding field evolution of the diffusion 
equation (^9|). Using the approximate evolution equation ( pl| ) together 
with a reduction operator 1Z of order M = 8 no difference occurs compared 
to the field evolution using all degrees of freedom. However for a reduction 
factor A = 2 as shown in figure very small deviations can be observed. 
We would like to point out, that both the exact calculation and the numer- 
ical results are highly dependent on linearity and adjointness of the Hamil- 
tonian, which account physically for the existence of normal modes. A 
similar calculation for non-self adjoint linear operators using singular value 
decomposition is given in the appendix. 
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Figure 2: The initial and the evolved fields for a lattice composed of 32 
sites. The approximate field evolution (□) is calculated according to the 
relation j2l| ) with a reduction of order M = 8 (a) and M — 16 (b) in 
the degrees of freedom. The field is plotted after 100 time steps using a 
temporal integration interval At = 1.0. The fields were evolved using a 
diffusion parameter v = 0.01. 

5 A General Approach to Nonlinear Evolu- 
tion Equations. 

In this section we apply the theoretical concepts developed in section || and 
section [|to the deterministic Kardar-Parisi-Zhang (KPZ) equation M and 
the deterministic Burgers equation |28| . Each equation provides a different 
type of non-linearity to the numerical minimization process and exhibits a 
different physical characteristics. As the main milestone in the direction 
of nonlinear surface growth the KPZ equation has been intensively inves- 
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tigated as the correct interface equation governing the physics of lateral 
growth phenomena p9[ . The Burgers equation describes a vorticity-free 
compressible fluid flow and the velocity field v(a:, t) becomes the analogue 
of the (surface) height function h(x,t) in the KPZ equation p0| |. 
In this section we monitor the evolution of both fields h(x,t) and v(x,t) 
over different time scales. We compute the scaling characteristics for the 
field fluctuations and various observables are measured in the long time 
scaling regime. We compare the performance of the evolution using a re- 
duced number of degrees of freedom with the evolution by direct numerical 
integration in time. In addition, by using transformations both equations 
can be mapped on the linear diffusion problem which was exactly solved 
regarding the reduction operator 1Z in section ^. These transformations 
provide an alternative exact construction of the reduction operator 1Z for 
both the KPZ and the Burgers equation. 

5.1 The relevant degrees of freedom in the KPZ equa- 
tion 

The deterministic KPZ equation is defined as 

= *,V 2 h(M) + 4Vh(M)] 2 , (41) 

where the first term on the right-hand side describes diffusive relaxation of 
the surface and the second term introduces the desired sideways growth. 

2G 



We discretize the linear part of equation ( [4ij ) as proposed in ( |30| ) and the 
additional non-linear evolution operator defined in (0) as 



«/(2At) if j = i-l;k = i-loTj = i + l;k = i+l 

:= \ -K,/(2At) 2 if j = i - 1; k = i + 1 or j = i + 1; k = i - 1 

else . 

(42) 

In analogy to figure [l] in section ^ figure || shows the result of the short- 



time RG approach introduced in section 3.1 as applied to the initial sur- 
face profile. The approximately evolved field displayed in figure || has been 
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Figure 3: The initial and the evolved surface profiles for a lattice composed 
of 32 sites after 100 time steps (a) and 300 time steps (b) using a temporal 
integration interval At — 1. The dashed curve represents the initial surface 
profile and the fields were evolved using a surface tension v = 0.01 and a 
growth velocity k = 0.1. 

evolved according to a reduction of order M = 16 in the degrees of free- 
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dom. Although the exact field evolution is recovered up to a time t = 100, 
for t = 300 the approximately evolved field exhibits a significant devia- 
tion from the evolution using all the degrees of freedom. For even larger 
times this results in uncontrollable numerical instabilities. The reported 
behaviour is expected for RG techniques which are not able to cope with 
scale interference pj| . 

To incorporate possible scale interference in the truncation procedure it is 



necessary to apply the concepts developed in section 3.2 and minimize \£\f 
in (|l^). In addition to the direct numerical minimization procedure we 
further propose a more sophisticated approach based on the exact solution 
previously established in section [I| There it was shown that for the linear 
evolution problem H = H the dynamic RSRG approach is equivalent to 
an exact diagonalization of the evolution operator. Since equation 
can be mapped to the linear diffusion equation (|2|) using the Hopf-Cole 
transformation j7| 

w(x,t) = e (.K/")M*,t) f (43) 

we are able to provide an exact reduction operator R to evolve the field 
under a reduced number of degrees of freedom []. In this case we decompose 
the construction of the reduction operator into three steps. First the initial 
field configuration h(a;,0) is transformed to w(x, 0) according to equation 

t Utilizing the Hopf-Cole transformation some authors reported about less instabilities 
in numerical integration schemes Jl6| 0|. 
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(p3|). The exactly calculated reduction operator including the relevant de- 
grees of freedom of the linear evolution operator H is applied within 
to evolve the field in time. Using the inverse of transformation (^3|) the 
final field is recovered by using a reduced number of degrees of freedom. 
In figure |J the final evolved growth profile h after time t — 100 is dis- 
played for lattices of different size. The approximate field evolution is 
performed using a reduction operator calculated by direct numerical mini- 
mization (square) of \£\f and alternatively using the exact transformation 
( (43| ) . Figure ||(a) shows the evolved growth profiles for a reduction factor 
of A = 4/3 for a periodic lattice composed of 16 sites, i.e. a truncation of 
order M — 4. In figure [|(b) the initial field is equally evolved using a reduc- 
tion factor A = 2. Compared to figure |^(a) the approximate equation 
fails to evolve the field correctly. In Figure [|(c) and (d) the correspond- 
ing surface evolutions are displayed for a lattice composed of 32 lattice 
sites. The results displayed in figure ^ demonstrate that the construction 
of the reduction operator 1Z is sensitive to finite size effects. In particular 
the reduction operator constructed by direct numerical minimization on a 
lattice composed of 16 sites fails to generate the characteristic paraboloid 
segments 0. 

Analogously to the numerical minimization procedure of the linear diffu- 
sion problem carried out in section ^, the SQP minimization algorithm 



was provided by the NAG Library 27 . Again the optimal solution was 
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Figure 4: The initial and the evolved fields for a lattice composed of 16 
sites ((a) and (b)) and 32 sites ((c) and (d)) after 100 time steps using a 
temporal integration interval At = 1.0. The dashed curve represents the 
initial profile according to a periodic random initialization. The fields were 
evolved using a surface tension v = 0.01 and a growth velocity k = 0.1. 
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found by the algorithm which in turn characterizes the minimization ap- 
proach in ( p7j ) as a well defined and numerically stable procedure also in 
the non-linear case where no general analytical approach exists. However 
the numerical minimization turns out to be much more time consuming 
since the error operator norm \£\f in equation is increasingly complex 
for non-linear evolution problems. 

To characterize the scaling behaviour of the KPZ equation we consider the 
decay of the density of surface steps for a lattice composed of N sites, 
defined by 

"(*) := ^E<i vh *wi> ■ ( 44 ) 

Here the brackets ( ) denote averaging with respect to an ensemble of initial 
surfaces (t = 0) characterized by the covariance 

(Ih.-h.l) ~ \%-j\< (45) 

with the roughness exponent £. The dynamic observable defined in equa- 
tion (|44|) obeys the scaling relation |Q 

p(t) ~ t ^- 1),z (46) 

where z denotes the dynamic exponent for the lateral correlation length 
£ determined by £ ~ t x l z . For the deterministic KPZ equation a scaling 
relation for z can be established in terms of £ as |3lJ 

z = 2 - C . (47) 
31 



To generate an ensemble of initial surfaces according to ( [45| ) we initialize 
the surface profile by the graph of a one dimensional Brownian bridge as vi- 
sualized in figure |^. The initial surface profile shown in figure ^| corresponds 
Graph of a periodic random walk Initial surface profile 



□ □ □ □ □ 
□ □ □ □ □ □ □ 
□ □ □ □ □ □ 

□ □ □ □ 



Number of steps 
of the walker 




Figure 5: An initial surface profile generated from the graph of a random 
walk using periodic boundary conditions (Brownian bridge). 

to a Brownian bridge of length TV introducing a finite size correction to the 
correlation function which vanishes for N — > oo. This dependence on 
the lattice size needs to be considered in a measurement for the roughness 
exponent of the theoretical value £ = 1/2 J35|. According to the techni- 
cal realization of the periodic boundary conditions the probability for the 
walker to go left and right after a certain number of surface steps depends 
on the remaining steps to go and the position of the walker. 
In figure |^(a)-(d) the density of surface steps observable is compared for 
four different reduction factors A = 4/3,8/5,2 and 8/3. We expect a de- 
viation from the exact scaling relation in the early time regime because 
of the finite size dependence of the initial surface profile constructed as a 
Brownian bridge. Therefore, inserting ( f47| ) into equation ( p6[ ) we derive 
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Figure 6: The inverse of the density of surface steps observable plotted 
according to the scaling relation ~ \ft for a different number of 

eliminated degrees of freedom (A = 4/3, 8/5, 2 and 8/3). The observable is 
measured over a total time of 100 time evolution steps of temporal distance 
At = 1 and the ensemble average was taken over 400 surfaces with £ = 1/2. 
The surfaces were evolved using a surface tension v = 0.01 and a growth 
velocity k = 0.1. 
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the scaling relation pit)' 1 ~ t 1 ^ 3 plotted in figure |] valid in the long time 
regime. As expected both curves show a significant deviation from the 
theoretical scaling behaviour in the short time regime. For longer times 
the evolution generated by using a reduction operator providing more than 
half of all the degrees of freedom for the nonlinear surface growth process 
recovers the desired long-time linear behaviour. Using a reduction factor 
of A = 2 or an even higher order in the reduction of the degrees of freedom 
the observable measured using the approximate field evolution displays an 
increasing deviation from the measurement using a field evolved with all 
the degrees of freedom. 

As a central observation in surface growth phenomena surface fluctuations 
exhibit a dynamical scaling behaviour J33| . The surface width of a growing 
surface is defined as [03l 



I N 2 I N 

\^E( h ^-^) with ^ = ]v^ h<,t - (48) 

\ i=l i=l 



In particular the scaling of the interface width ( |4q ) on the length scale N 
is expected to be of the form |29| 

W t = Ntf t/N . , (49) 

where ft/N z denotes a scaling function different from the one in the stochas- 
tic case. Using relation ( ff7| ) we calculate for the KPZ equation the scaling 
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characteristics 

1 



N m = ft/m/i ■ (50) 



In figure [?] the surface width observable is plotted in analogy to figure || 
measured using equally evolved surface profiles. All measurements of the 
surface width observable are calculated according to an ensemble average 
of 400 surface configurations to be consistent with the literature [[H] . The 
measurements for the observable W represent the same dependence on the 
reduction factor A as the analogue measurements for the observable p in 
figure ^[ 

5.2 The relevant degrees of freedom in the Burgers 
equation 

As a quasi one dimensional analogue of the Navier-Stokes equations the 
Burgers equation is defined as 

gV ^' ^ = vV 2 v(x,t) ~ Kv(x,t) ■ Vv(x,t) , (51) 

where the parameter v denotes the viscosity of the velocity field v(x,i). 
The nonlinear term on the right hand side is a transport or convection 
term in which the speed of the convection depends on the magnitude of 
v(x, t) 0. Compared to the KPZ equation analyzed in the previous section 

tin the original equation as introduced by Burgers k = 1 Pq . 
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Figure 7: The surface width observable plotted according to the scaling 
characteristics of the KPZ equation. From (a)-(d) increasing reduction 
factors of A = 4/3, A = 8/5, A = 2 and A = 8/3 are used. The observable is 
measured over a total time of 100 time evolution steps of temporal distance 
At = 1. The ensemble average was taken over 400 surfaces with £ = 1/2 
each generated by a periodic random walk as displayed in figure |^. All 
surfaces were evolved using a surface tension v = 0.01 and a growth velocity 

K = 0.1. 
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the convection term provides a different discretization scheme to the direct 
numerical minimization procedure. 

Inserting the transformation v(x, t) = (d/dx)h(x,t) into equation ( |5l| ) we 
rederive the KPZ equation defined in (f4l|). By further applying (^3|) results 
in the linear diffusion equation for which we can calculate the reduction 
operator exactly as described in section ||. Again, following the steps out- 
lined for the case of the KPZ equation below (|43| ) we are able to provide an 
alternative exact reduction operator 1Z to evolve the velocity field under a 
reduced number of degrees of freedom. 

Figure || shows an evolved velocity field after a total evolution time i = 100 
using a lattice composed of 32 sites. In figure ^ the performance between 
the reduction operator 1Z constructed by directly minimizing \£\f and the 
exact one using the proposed transformations is compared for two different 
reduction factors A = 4/3 and A = 2. Although both methods perform rea- 
sonably well for a reduction factor A = 4/3 for a reduction factor A = 1/2 
the direct numerical minimization performs far superior to the evolution 
incorporating the exact transformations. The stated performance of the 
SQP algorithm using the NAG Library f27|] is equal to the examined case 



of the KPZ equation in section 5.1, so that also in this case the numeri- 



cal minimization converges in a well defined minimum subject to equation 



Burgers showed that the decay of the step density follows the dynamical 
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Figure 8: The initial and the evolved velocity fields for a lattice composed 
of 32 sites after 100 time steps using a temporal integration interval At = 1. 
The initial velocity field configuration (dashed curve) is evolved with a re- 
duction factor A = 4/3 (a) and A = 2 (b). The evolution using a reduced 
number of degrees of freedom is shown using a reduction operator con- 
structed by numerical minimization and by using an exact transformation. 
The velocity fields were evolved with a surface tension v = 0.01 and a 
growth velocity k = 0.1. 



scaling behaviour [p8| 



<|v(M)l) ~ ts 



(52) 



which is equivalent to the scaling relation (fi"6| ) with £ = 1/2 and z = 3/2 
using the transformation v(x,t) = (d/dx)h(x,t). Here the brackets () 
denote an ensemble average of stationary and locally uncorrelated initial 
velocity field configurations v(x, 0). In figure |^ we have plotted the asymp- 
totic time behaviour over a total time t — 1000. The velocity field 
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Figure 9: The inverse of the decay of the density of surface steps observable 
plotted according to the scaling relation ( |52|) for a reduction factor A = 4/3 
(a) and A = 2 (b). The observable is measured over a total time of 1000 
time evolution steps of temporal distance At = 1 . The ensemble average 
was taken over 400 random surfaces. All surfaces were evolved using a 
surface tension v — 0.01 and a growth velocity k = 0.1. 

has been evolved by a reduction of order M — 8 and M = 16 degrees of 
freedom respectively, the latter corresponding to a A = 2 reduction. Fol- 
lowing Burgers j28| we generated an ensemble of stationary and locally 
uncorrelated initial velocity field configurations Vi(a;,0) by calculating 



Vi(x,0) = (d/dx)hi(x,0) and % e {1, . . . ,400} 



(53) 



In equation (|53] ) hi(x, 0) denotes the graph of a Brownian bridge as shown 
in figure ||. 

To confirm the results in figure || and ^ quantitatively we investigate the 
numerical stability within the evolution of the velocity field in dependence 
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of the order of the truncation in the limit of long-time evolutionary dy- 
namics. Using ( pl| ) together with an iterative construction of the reduction 
operator 1Z M of order M according to equation (|2^) , table [I] summarizes 
the mean square deviation (MSD) as measured in comparison with the ex- 
act velocity field evolution. The convergence in time of the MSD for each 

Table 1: The accuracy of the approximate velocity field evolution using 
direct numerical minimization as compared to the exact velocity field evo- 
lution using all degrees of freedom in the long-time limit. The approximate 
velocity field evolution is compared for a reduction of different orders by 
calculating the mean square deviation (MSD) between the approximately 
and the exactly evolved velocity field. All reductions calculated refer to an 
original lattice composed of 32 sites. 
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order of the reduction in the degrees of freedom illustrates the numerical 
stability in the long-time scaling regime. However, for a reduction factor 
A > 2 (M = 16) the accuracy of the approximately evolved velocity field in 
the long-time evolution decreases by at least two orders of magnitude. This 
indicates that at least half of the available degrees of freedom are relevant 
for an accurate evaluation of the long-time characteristics in the Burgers 
equation. 

6 Conclusions 

In this paper we applied an operator real-space RG method to partial 
differential equations (PDEs) in general. The fundamental concept is to 
construct a linear map, the reduction operator 1Z, projecting to a subspace 
of the full vector space including the relevant degrees of freedom governing 
the physics of the evolutionary dynamics. In technical terms the reduction 
operator is composed as the product of the truncation and embedding op- 
erator, which are already established objects in the equilibrium real-space 
RG literature. 

Two basic approaches are provided for the construction of the reduction 
operator 1Z. In the first approach a transformation to a linear PDE is used 
for which the reduction operator can be constructed using exact or numer- 
ically fast diagonalization techniques. If such a transformation is missing 
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a direct numerical minimization is applied using a sequential quadratic 
programming (SQP) algorithm to construct the reduction operator. The 
advantage of the first approach over the second one is dominated by the 
saving of computer time depending on the particular nonlinear problem of 
interest. However, again depended on the particular PDE, it was observed 
in this work that numerical instabilities in the field evolution can prove 
the first approach to be useless. Although the computer calculation time 
increases with the lattice size the computed reduction operator does not 
depend on a particular initial field configuration and can be stored and 
applied to the evolution of any other initial field configuration. 
In fully developed turbulence or related non-integrable dynamical systems 
the dynamics and evolution of many degrees of freedom interacts through 
nonlinear PDEs. Direct numerical simulations are difficult and show a clear 
demand for a reduction to the relevant degrees of freedom. Choosing the 
Burgers equation to test our RG approach we have examined one of the 
simplest archetypes describing a vorticity-free, compressible fluid flow. The 
measured decay of the step density can be generalized to higher order cor- 
relation functions corresponding to structure functions in fully developed 
turbulence fl30f . 

One of the goals of the real-space RG method presented in this paper is 
certainly the reduction in complexity of the problem yielding to a reduced 
amount of computer time within simulations, after a reduction operator 1Z 
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has been provided. Since this operator contains information on the rele- 
vant degrees of freedom the proposed real-space RG method gives rise to 
improved understanding of the physics of non-linear evolutionary dynam- 
ics. 
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Appendix 

In the case of a non-selfadjoint Hamiltonian, the approach in section ^| must 
be slightly changed. There exists no decomposition of the functional space 
into eigenvectors of the evolution operator H , so the natural approach re- 
quires a singular value decomposition pj| 

JV 

H = J2 E i\ u i)( y i\ ■ ( 54 ) 
<=1 



In equation (54) the vectors denote the "input states" and the vec- 

tors make up the "output states" . It is important to state that both 

sets of vectors belong to the same functional space. 
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Changing slightly our notation for clarity \<j>) denotes the target state, 
whose removal minimizes the error operator 

£ := P W + At [HV W + V W H - (0| H |0) P w ] . (55) 

We decompose \cf>) in both the input and the output basis as 

N N 

\4>) = ^2^\v t ) and \<j>) = ^XiK) . (56) 

i=l i=l 

Analogously to (Q) we arrange the resulting operators in (|55| ) using the 
two alternative representations in ( |56| ) by 

N N 

V \4>) = X XtMj \ u i)( v j\ (<t>\ H \<t>) = Y^iXiEi = (H) 

i—l 

N N 

HV \4>) = X] E i^o K) ( v i\ V \4>) H = X] XiXjEj\ui){vj\ . 



(57) 



Therefore, in the basis of the evolution operator 

N 



£ = £ij \ U i) ^ « 



with £y = XiMj + At^Eifiifij + XiXj E j ~ ( H ) Xt^j) ■ (58) 
Squaring and retaining only terms which are linear in At 

4 « X*/*? + 2At(E ilHXi $ + EjxlxjH - (H) X^f) • ( 59 ) 

Summation over all the values in |[|) yields 

iv 

£ 4 = 1 + 2-At{H) (60) 
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to first order in At. The result is equivalent to equation ( |37| ) and a again 
(H) must be minimized. However one has to remember that the set of 
numbers {/ii}; and {xi}i are not independent since they represent the same 
field. 

The idea is to represent the vectors in the basis of the vectors 

according to 

N 

The matrix CV,- may be obtained by the SVD of the evolution operator H 
using 

Cij = \vj) ( Ui \ . (62) 



The equality of both representations in ( |56| ) yields a relation between the 
and the {xi}i- 



N N N / N \ JV 

X)Mi|«i> = Yl K) = [^2 Ci ^ij \ u i) = ^^luj) (63) 

i— 1 id = 1 J = l \i=l / J — 1 

As the {|itj)}i form an orthonormal basis we may deduce 

N 

Xj = ^djVi ■ (64) 

i=l 

Using (^) and inserting ( |6~4| ) the value to be minimized equals 

JV N 

(H) = (0| if |0) = ^mxiEi = J2 E ^i°3iN ■ (65) 
4—1 i 
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Equation ( pq ) defines a quadratic form based on a matrix which is not 
symmetric. Defining 

Kii^^EiCji + EjCa) (66) 
the problem is reduced to the minimization of the quadratic form 

N 

(H) = ]T KiiMj (67) 

with the restriction ^ /if = 1 . This problem is equivalent to the self- 
adjoint case examined in section ^ and the solution is the lowest eigenvalue 
of the matrix K. 

The eigenvector corresponding to the minimum eigenvalue of the matrix K 
is to be interpreted as a set of values. Thus, the \<j>) eigenvector shall be 
given by the calculation of J2i Mi \ u i) with the appropriate values of 
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